initialization

suppressMessages({
  ddPath <- "/Users/michaelalfaro/Dropbox/color_evolution_scratch/tanager_comp_analysis/allison_data/RefDataSpeciesSummary_MatchingMF"

  dat <- pblapply(list.files(ddPath, full=T), read.csv)
  wl <- dat[[1]][, 1]
  dat <- sapply(dat, "[", -1)
  
  specs <- as.rspec(cbind(wl, do.call(cbind, dat)), lim=c(300, 700))
  specs <- procspec(specs, fixneg="zero")
  
  male_specs <- specs %>% select(!ends_with("f"))
  
  crown_data <- male_specs %>%
    select(wl, starts_with("Crown_")) %>%
    rename_with(~ gsub("m$", "", .)) %>%
    rename_with(~ gsub("^Crown_", "", .))
  
  rgb_colors <- spec2rgb(crown_data)
  names(rgb_colors) <- names(crown_data)[-1]
  
  head(crown_data)
})
##    wl   AcaBai    AniIgn   AniLac    AniMel   AniNot   AniSom   BanArc   BanAur
## 1 300 2.356690 1.1447721 11.39746  8.896685 4.535676 2.014334 5.533874 3.126775
## 2 301 1.903805 0.8218703 12.11596  8.899757 5.386252 1.444982 4.656108 3.745577
## 3 302 1.777677 0.3635018 12.44752  9.508270 6.579554 1.646972 4.473581 4.624964
## 4 303 1.956667 0.5506667 12.46547  9.872667 6.587833 1.983947 4.986000 4.196333
## 5 304 2.086104 0.5410315 12.57834  9.691481 6.914120 2.683912 5.666046 5.145111
## 6 305 2.553706 0.9282766 12.29238 10.253288 7.446329 2.923659 5.490410 6.381721
##      BanEdw   BanMel   BanRot   ButAur    ButExi ButMon   ButWet CalCoc
## 1 1.4074342 2.792863 2.720081 12.36975  9.853977      0 1.940090      0
## 2 0.8987856 2.613824 2.721775 12.86508  9.873162      0 2.141559      0
## 3 1.0957351 2.701293 2.913748 13.68989 10.111072      0 2.995856      0
## 4 1.6188667 2.551333 2.243000 13.83140 10.321933      0 2.810667      0
## 5 2.0386333 2.740259 1.734444 13.45717 10.115626      0 3.157426      0
## 6 2.5069297 3.211658 1.078811 14.04848 10.483344      0 3.042622      0
##     CamHel   CamPal   CamPar   CamPau   CamPsi   CatAna   CatDia   CatHom
## 1 3.469025 6.643459 2.719078 5.499726 2.579595 4.518847 2.219834 5.587705
## 2 3.036942 7.602843 3.744938 5.698823 3.010506 4.336070 2.382877 6.276141
## 3 2.725917 7.989854 4.219054 5.927595 3.078775 4.244808 2.497530 6.854827
## 4 2.889133 7.241067 3.745280 5.383200 3.367667 3.862467 1.883267 6.427067
## 5 3.143063 6.617119 3.128248 5.333707 3.731115 4.145883 2.269567 5.481470
## 6 3.869150 6.013214 3.264869 5.477861 3.516852 4.021309 1.918708 6.225052
##     CatIno   CerOli   ChaEuc   ChlCal ChlNit   ChlPho   ChlRie   ChlSpi
## 1 3.008490 1.597151 5.387928 2.630297      0 1.787583 2.784198 1.612941
## 2 2.616133 1.730229 6.023412 3.474346      0 2.044171 3.276139 1.277324
## 3 3.294292 2.450377 6.458176 3.473768      0 1.696432 3.689607 1.380377
## 4 4.443133 2.863067 6.054750 3.991000      0 2.461333 4.117933 1.372000
## 5 3.437819 2.587056 5.016699 3.586952      0 3.661704 4.370585 1.306522
## 6 5.182589 2.042841 5.806554 3.673250      0 3.939279 4.340295 1.892305
##     ChrChr   ChrSal   CisLev   CneRub   CoeFla   ComBae   ComGar   ComLor
## 1 1.949883 4.471500 6.664564 4.975148 3.760420 4.967959 1.844345 4.560399
## 2 2.329393 4.400425 6.275407 4.706490 4.126041 5.326747 2.188339 4.799556
## 3 2.523396 4.784778 6.921449 5.457355 4.167218 4.105865 1.724850 5.719318
## 4 2.792556 4.523722 6.041067 6.047067 4.111000 3.140896 2.095556 5.047583
## 5 2.816093 3.580108 6.493004 6.208793 4.068874 3.156895 2.492580 4.814028
## 6 2.818090 3.429653 6.855703 6.840187 4.819775 2.791245 2.552961 4.056009
##     ConAlb   ConBic   ConCin    ConFer   ConLeu   ConMar   ConRuf   ConSit
## 1 10.46993 2.015616 3.804209 1.9168577 2.578745 4.845310 3.141086 2.525085
## 2 10.12310 1.258997 3.714779 1.8302739 2.849072 4.560984 2.548177 2.491962
## 3 10.40977 1.210195 3.284392 0.9475532 3.220024 4.541942 2.517834 2.725701
## 4 11.18733 1.583111 3.498333 0.9898000 2.851333 4.683467 2.070467 2.394800
## 5 11.26497 1.388296 4.440231 1.8609481 2.321988 4.088907 2.744841 2.564296
## 6 12.11861 1.712550 4.598439 1.9402180 1.952958 3.663535 2.913573 2.412890
##     ConSpa ConSpm   ConTam   CorCuc   CorMel   CorPil   CreDen   CreVer
## 1 3.358489      0 3.412506 6.490117 4.495232 10.35447 2.840637 4.538499
## 2 3.709069      0 3.201041 7.144863 4.478360 10.73900 2.686527 4.503564
## 3 3.925063      0 3.037099 7.932816 4.588965 10.93741 1.949874 4.403225
## 4 4.112333      0 3.636267 7.416467 4.771722 12.02407 1.951083 5.057800
## 5 4.758173      0 4.172085 7.681907 5.639179 10.85867 1.922611 4.922030
## 6 4.390084      0 4.023079 7.862213 5.623269 10.67903 1.630532 5.197043
##     CyaCae   CyaCys   CyaLuc   CyaNit    CyiCyi   CypHir   DacAlb   DacBer
## 1 47.82842 41.49798 23.52521 36.24909  9.682599 3.422173 30.74110 10.52712
## 2 48.71184 42.43106 24.20018 37.04390 10.122437 3.426750 30.94648 11.42346
## 3 49.29641 43.49012 24.71517 37.55738 10.289468 3.745732 31.15261 10.78989
## 4 49.57293 44.47547 24.95913 37.89383 10.688333 3.887600 31.72800 11.02356
## 5 50.42770 45.24845 25.83192 39.48137 10.758269 3.550522 32.22233 11.77883
## 6 50.60203 45.97401 26.38695 40.60357 10.539698 3.741910 32.99691 11.63292
##     DacCay    DacFla   DacLin   DacNig   DacVen   DacVig   DelCas   DigAlb
## 1 16.62614  7.034602 22.20075 20.04237 24.05480 20.14063 12.71208 5.808184
## 2 17.32347  7.183602 22.53657 20.95381 24.02077 20.86264 13.64604 5.711144
## 3 17.80806  8.018494 23.07261 20.84107 23.55332 21.15301 13.70460 5.501305
## 4 18.50593  8.506333 24.06907 20.39833 23.70200 21.38573 13.86287 5.799800
## 5 19.93100  9.992485 24.85467 20.41631 24.27682 21.83977 14.96010 6.483400
## 6 19.51759 10.959667 26.11854 21.03806 24.74464 22.36942 15.42640 6.707083
##     DigBar   DigBru   DigCae   DigCar   DigCya   DigDui   DigGla     DigGlo
## 1 6.442543 5.176213 4.761571 3.826953 10.33006 7.791004 12.56828 0.00000000
## 2 5.508182 5.402843 5.634023 3.892699 10.60421 8.172086 12.80130 0.00000000
## 3 5.339232 5.436029 5.975823 2.957868 10.81254 7.808751 12.78159 0.03230631
## 4 4.667250 5.552800 6.400200 3.486667 11.47007 6.575733 13.07420 0.00000000
## 5 5.529088 5.609878 6.701107 4.228474 11.73469 6.108693 13.32874 0.41992963
## 6 5.729264 5.660191 6.361795 3.445634 12.32895 5.910584 13.70849 0.95878018
##     DigGls   DigHum   DigInd   DigLaf   DigMaj   DigMys    DigPlu   DigSit
## 1 5.911885 5.315526 3.912009 4.740884 5.746431 6.459852  9.968993 7.514505
## 2 5.905845 5.477690 3.296104 5.007212 5.998514 5.987995 10.236485 8.075579
## 3 6.023903 5.611667 3.129428 4.902143 6.072998 6.033641 11.279784 7.875029
## 4 6.213250 5.820533 2.860333 5.006771 6.509533 6.716333 11.134667 8.076500
## 5 6.421384 6.156630 2.969861 4.906539 5.850556 6.852841 11.711237 7.907699
## 6 6.629507 6.380209 3.943491 4.563189 5.928620 7.607613 12.199932 7.886867
##     DigVen   DiuDiu   DiuSpe   DolFri   DonAlb    DubTae   EmbDui   EmbHer
## 1 4.413730 5.754768 4.513088 5.628032 4.645562  9.388712 1.686559 4.040586
## 2 4.292516 6.051697 4.685247 6.604378 5.038362  9.714932 2.270955 3.951989
## 3 4.044523 5.380742 4.342351 7.106216 4.673121 10.139378 1.531207 3.911493
## 4 3.733833 4.769333 4.060133 6.991333 4.297933 10.404267 3.147333 4.198333
## 5 4.270157 3.463952 4.665000 6.549444 3.256133 10.792793 3.227259 3.965176
## 6 4.472887 2.724789 4.475528 5.723248 2.586216 10.657593 3.542144 4.317468
##     EmbLon   EmbPla   EmbYpi    EucPen    EunCam   GeoCon   GeoDif   GeoFor
## 1 2.338532 4.269550 4.909153 0.4650937  8.777310 4.302198 4.160865 5.656685
## 2 3.439315 3.636867 4.307609 0.8571243  9.311593 3.814930 4.490456 5.751912
## 3 3.362640 4.218802 3.409551 1.2577694 10.058741 4.273697 4.224193 6.028177
## 4 3.698667 3.939067 4.183800 1.4741333 10.300400 3.505800 4.474000 6.035600
## 5 3.461648 3.932163 5.405856 1.1612111 11.166744 3.235274 4.708181 6.405648
## 6 3.409676 4.654184 6.260063 0.9784072 10.824191 3.115571 4.505560 6.617688
##     GeoFul   GeoMag   GeoSca   GubCri   HapRus   HapUni   HemAtr   HemCal
## 1 6.448760 6.531546 4.561547 7.380761 4.582564 7.463453 2.444315 1.789450
## 2 5.698004 7.064521 4.258554 6.909324 5.085793 7.039780 2.618054 1.943099
## 3 5.313429 7.979640 4.555024 6.311311 5.553897 8.032144 2.213946 2.463144
## 4 5.100000 8.089267 3.589333 6.047167 5.761733 8.117767 2.134833 2.643067
## 5 4.763537 8.176015 3.280883 5.714065 5.783456 8.154983 2.177606 3.226993
## 6 5.055092 7.564323 2.848854 6.350489 5.951241 7.519610 1.680514 3.041022
##     HemFla HemFro   HemGoe   HemGui   HemMel    HemPar HemPiu   HemRey   HemRuc
## 1 3.287751      0 5.605915 3.379796 3.736193 3.1149820      0 1.266595 4.524636
## 2 3.218171      0 5.699159 3.965413 3.895222 3.4821216      0 1.997243 4.850843
## 3 2.946195      0 4.831592 4.463416 3.758582 1.4672354      0 2.621126 5.460076
## 4 2.740067      0 4.679933 4.512533 3.976467 0.9297917      0 3.241667 5.154133
## 5 2.694589      0 4.609824 4.911848 4.363937 0.8993449      0 3.903537 3.785207
## 6 2.550214      0 4.942499 3.846759 4.760578 1.1250101      0 4.106505 3.929441
##     HemRus   HemSup   HemTri   HemVer   HemXan   HetRub   HetXan   IdiBra
## 1 2.502548 2.098243 3.904831 3.234018 3.045542 6.152712 4.228811 3.312871
## 2 2.309180 1.753862 4.316602 3.620662 2.942859 6.544171 4.682160 3.618069
## 3 2.877843 2.183583 5.085332 4.342088 3.038568 6.654640 4.674836 4.216775
## 4 3.275933 2.895000 5.113800 4.343833 3.104200 6.115667 4.887733 4.643333
## 5 3.262852 2.984352 5.068830 4.830097 3.556463 6.642037 4.921104 4.764395
## 6 3.059625 3.220811 5.025465 5.040315 3.916326 6.796928 4.661458 4.698480
##     IncLae   IncOrt   IncPer   IncPul   IncWat   IriAna   IriJel   IriPor
## 1 4.414038 5.292791 2.551869 3.270497 4.785303 5.771081 3.992751 4.377826
## 2 3.989401 6.119745 2.390995 3.248391 5.037771 5.888508 4.252169 3.813357
## 3 4.258971 6.272973 2.418378 3.395009 5.340468 6.093108 4.324804 4.009715
## 4 3.983333 6.503333 2.145333 3.724733 5.192800 5.781933 4.280600 4.185222
## 5 3.982542 6.122935 1.942009 4.036441 5.664348 5.276085 4.193115 4.397148
## 6 4.378187 5.115919 2.175338 4.137515 5.404914 4.974350 4.235859 4.410030
##     IriPul    IriRei   IriRuf   LanAur   LanFul    LanLeu   LanVer   LopGri
## 1 2.556881 1.6641640 2.398438 1.711784 2.155905  6.806644 1.875807 4.396759
## 2 2.662247 1.6932523 2.305449 2.010827 2.309213  8.114988 1.914119 4.427661
## 3 2.329413 1.1432036 2.183699 1.496450 2.127450  8.692264 2.270701 4.317014
## 4 2.325067 1.1347333 2.418600 1.369267 1.809733  9.683662 2.384200 3.895667
## 5 2.009063 1.2780593 3.164674 1.538796 1.728348 10.244642 2.234641 3.910170
## 6 1.805400 0.9912414 3.542218 1.353515 1.455249  9.401372 2.153353 4.449196
##     LopPus   LoxAno   LoxNoc    LoxPor   LoxVio   MelMel   MelNig   MelRic
## 1 7.488234 3.910953 3.112941 0.7241694 3.978859 6.290421 3.317616 2.663009
## 2 7.728872 3.611694 3.302369 0.5207532 3.494323 6.346495 3.616895 2.822505
## 3 8.120056 3.903131 3.412677 0.4182378 2.825022 5.101267 3.853908 2.791955
## 4 8.777200 3.814917 3.506333 0.8708667 2.716667 4.902100 3.827133 1.614333
## 5 9.125893 3.826162 3.245481 0.7432741 2.436115 4.756217 4.466144 2.336722
## 6 9.077670 4.258948 2.971335 0.7415405 2.341905 4.632594 4.545243 2.238441
##     MelXan   NemPil   NeoFas   NepOne   NesAcu   OrcAbe   OreFra   OryAng
## 1 3.227230 4.226787 5.055903 2.750476 4.476550 4.686505 4.003340 4.710562
## 2 3.053802 4.427649 5.578859 2.984126 4.927500 5.367763 4.011671 4.578402
## 3 3.187752 5.991195 5.745822 3.277987 5.451937 6.685174 3.800905 4.211364
## 4 3.999750 6.571733 5.922933 4.008067 5.517333 5.427667 4.056417 3.453000
## 5 4.214352 6.875541 6.054052 3.772852 5.130167 5.959997 4.003097 3.298693
## 6 3.863957 7.280532 5.968306 3.802959 5.445279 5.704658 3.699477 3.357095
##     OryCra   OryFun   OryMax   OryNut   ParCap   ParCor   ParDom   ParGul
## 1 7.429500 3.171680 2.478050 3.547213 3.131365 2.192040 3.488714 5.262649
## 2 7.114906 1.696083 2.238608 3.181813 3.250432 1.692593 3.211316 5.341650
## 3 6.791849 1.783286 2.637338 3.065804 2.936459 1.480622 2.503977 5.598689
## 4 6.554433 2.014250 2.274000 3.493822 2.827667 1.769867 2.564167 5.632500
## 5 6.204209 1.687894 2.008000 2.797590 2.573741 2.164226 3.877600 5.221044
## 6 6.324417 2.340734 2.263856 3.435716 1.892505 2.509106 4.650066 5.927428
##     ParHum   ParNig    PhrAla   PhrAtr   PhrCar   PhrDor   PhrEry   PhrFru
## 1 1.291459 2.886041  9.293345 3.329982 4.506378 5.958768 5.393117 4.035439
## 2 1.470869 3.250788  9.722730 3.493220 4.852748 6.540127 5.173101 3.700482
## 3 2.591761 3.799027  9.971624 3.316380 5.439523 7.434704 5.215481 3.711192
## 4 3.333333 3.780667 10.363583 3.523467 5.290667 6.914625 5.228400 4.010433
## 5 1.288852 4.100083 11.004463 3.808385 6.077074 7.631400 5.338504 4.176483
## 6 2.063306 3.766527 11.143282 4.008009 5.946428 7.364779 5.950004 4.193376
##     PhrPat   PhrPle   PhrPun   PhrUni   PieCin   PinIno   PipMel   PlaCra
## 1 5.918407 4.492483 4.337694 4.419768 4.849804 7.595793 44.03170 4.935750
## 2 4.687977 4.608919 4.675339 4.146243 5.063029 7.981541 45.27175 4.996324
## 3 3.553551 4.704658 4.880333 4.382067 5.183631 6.137357 46.54217 5.160751
## 4 3.805467 4.587800 4.570267 4.332867 5.437083 5.480333 46.97613 5.289200
## 5 4.819904 4.487626 4.711126 3.984511 5.364856 4.238993 47.44355 5.450237
## 6 5.570063 4.528342 4.649456 4.118596 5.041036 3.369950 47.98022 5.174789
##     PooAlt   PooBol   PooCab   PooCae   PooCin   PooEry   PooHis   PooHyp
## 1 2.659149 5.074505 3.147468 1.979581 3.787056 7.712095 2.747193 4.714058
## 2 3.038234 4.900971 3.120582 2.135532 4.244384 7.609276 2.777488 4.660856
## 3 3.145950 4.980750 3.085348 2.241360 4.363910 7.682692 2.762346 4.884816
## 4 3.486833 5.398733 3.029067 2.709000 4.223600 7.945467 2.773867 5.114267
## 5 3.895370 5.529215 3.251222 3.093833 4.675441 8.458844 3.236433 4.961119
## 6 3.539613 5.523874 3.204924 2.964775 4.821710 8.432375 3.034631 4.787085
##     PooMel   PooNig   PooOrn   PooTho   PooTor   PooWhi   PorCae   PyrRuf
## 1 11.26016 5.275041 5.698532 4.047944 1.994281 6.175694 8.884247 2.985196
## 2 11.80131 5.415951 5.689101 4.329250 2.436692 6.643366 9.345823 2.704384
## 3 11.56038 5.513935 4.807838 4.291748 3.602207 6.608038 9.545780 3.434103
## 4 11.83930 4.856200 3.447583 4.118200 3.506933 6.874067 9.865067 3.132467
## 5 10.55111 3.877707 4.733398 4.123044 3.973230 6.915759 9.460619 4.185693
## 6  9.37007 4.818914 4.960957 4.444984 3.475429 6.550272 9.483804 3.727236
##      RamBre    RamCar     RamCos    RamDim    RamFla   RamMel    RamNig
## 1 3.8178378 0.6575369 0.00000000 0.0000000 0.0000000 2.373670 0.8585081
## 2 1.7594474 1.2629982 0.00000000 0.4934586 0.0000000 2.212339 0.7717676
## 3 2.2941652 1.8017676 0.03942703 0.9421703 0.1915694 1.945461 1.1554324
## 4 0.7708889 2.4197333 0.28406667 1.0868333 0.0000000 1.657600 1.4081333
## 5 0.5567407 2.7865519 0.34956667 0.8162907 0.2692222 1.029870 1.4319370
## 6 2.3220991 2.3640991 0.31947207 0.3772622 0.0000000 1.078897 1.7661766
##     RamPas   RamSan   RhoCru   SalAlb     SalAte   SalAto   SalAtp   SalAur
## 1 2.092480 1.107279 5.267655 3.915486 0.00000000 4.418674 3.056477 4.420184
## 2 1.974804 1.321750 5.313646 3.759941 0.00000000 4.314685 2.939376 4.057112
## 3 1.886745 1.813571 5.048610 4.188431 0.00000000 4.326058 2.848426 4.018236
## 4 1.575917 2.131267 4.617833 3.978600 0.00000000 3.803733 2.996133 3.447733
## 5 1.236782 2.613637 4.703222 3.953630 0.04198148 3.729526 3.649130 3.590822
## 6 1.318827 2.669377 4.916980 4.746939 0.53640541 3.898820 3.974416 4.002832
##     SalCin   SalCoe    SalFul   SalGro   SalMal   SalMam   SalMul   SalNig
## 1 3.168025 4.004351 1.0794342 5.033497 5.002832 5.334402 15.74235 3.370532
## 2 2.839908 3.822681 1.0258180 4.819243 4.476708 5.543121 16.62963 3.105339
## 3 3.011948 3.490829 0.5417351 4.965391 4.330586 5.199097 15.40084 3.240399
## 4 3.607250 3.170400 2.0172667 4.799333 4.866000 4.940267 15.30760 3.271667
## 5 3.169894 3.324511 2.5264296 4.767463 6.317667 5.091256 14.95960 3.644216
## 6 3.618851 3.915580 2.5987027 4.726150 6.961162 5.189872 14.91429 4.024294
##     SalOre   SalRuf   SalSim   SalStr   SchMel   SchRuf   SerAlb   SicAur
## 1 5.463356 5.274351 4.474921 6.131899 4.324312 4.449703 14.00662 8.577032
## 2 7.085903 5.478189 4.156344 6.256351 4.736539 5.146381 14.34296 8.028575
## 3 6.254764 5.436554 4.238843 6.078312 4.933793 5.448093 14.93730 7.989236
## 4 5.441367 4.852167 4.596400 6.137133 4.674733 5.316222 15.42970 6.809067
## 5 5.154061 5.456111 4.878281 6.368189 5.512852 5.474654 16.32746 5.878826
## 6 4.032126 5.520279 5.259483 6.556189 5.823301 5.291928 16.89555 6.798753
##     SicCit   SicCol   SicFla   SicLeb   SicLua    SicLuc   SicLuo   SicLuv
## 1 11.25139 7.205077 2.830766 7.648045 4.676791 0.1021811 5.362960 5.264173
## 2 12.29999 7.287677 3.230378 7.145766 5.741486 0.2458739 5.011879 5.571210
## 3 12.99739 7.838106 4.078973 8.224568 6.527688 0.6680757 4.439155 4.883621
## 4 13.91120 6.471133 4.317333 7.508000 6.496067 1.2284667 4.784400 3.816833
## 5 14.09174 7.217593 4.263859 7.562741 5.551548 1.5280704 5.216737 3.680078
## 6 14.09012 6.273910 4.061829 7.912486 5.579627 1.6636676 5.659854 3.144914
##     SicOli   SicRai   SicTac   SicUro   SpoAlb   SpoAme   SpoArd    SpoBoe
## 1 5.204306 6.755234 4.564286 3.310614 4.268679 2.761505 3.879212 0.2185383
## 2 5.255086 6.567703 4.546519 3.668602 4.282187 3.020553 4.298396 1.2393041
## 3 5.441187 5.799820 4.824760 4.187177 4.633883 3.492620 4.489554 1.7047027
## 4 5.620867 6.934000 4.833533 4.504667 4.380400 3.838067 3.044667 1.8463333
## 5 5.716030 8.905722 5.345681 4.650185 4.569707 2.745678 2.185917 1.6927407
## 6 6.155941 8.194748 5.210966 4.254263 4.497658 2.635780 1.390338 1.2230000
##     SpoBoo    SpoCae   SpoCas   SpoCin   SpoCol   SpoCor   SpoFro    SpoHyc
## 1 4.993282 1.9897946 4.061829 6.635216 4.723308 4.821036 4.289338 10.106671
## 2 4.729613 1.4693568 3.707445 7.844453 4.328332 5.098705 3.536815  9.856275
## 3 4.072123 1.2005856 3.342933 8.342039 4.142750 7.214793 2.711018 10.156014
## 4 4.578778 0.9673333 3.272533 8.349444 4.217000 6.789800 3.588333 10.500500
## 5 5.105025 1.6696407 3.542422 8.371451 3.992930 7.721230 2.995056 10.530398
## 6 5.205201 1.9081766 3.553863 8.365844 4.160517 6.525811 3.909784 10.799360
##     SpoHyx    SpoIna   SpoInt    SpoLeu   SpoLin   SpoLuc   SpoMin   SpoNic
## 1 9.887982 11.759180 4.064547 13.144916 5.005912 4.448468 4.661175 5.567274
## 2 9.209042 10.563775 4.201367 12.932435 4.390182 4.509087 5.045897 5.589960
## 3 9.126309  9.028432 3.419466 12.430396 3.883852 4.325381 5.108537 5.457667
## 4 9.069533  9.184333 3.452333 10.863333 4.033333 4.926000 5.080867 5.027067
## 5 7.520901  9.913889 3.455208  9.510302 3.316396 4.462105 5.410693 5.467456
## 6 8.636017 11.347477 3.394881  9.717895 3.347627 4.186982 5.259308 6.098886
##     SpoPer   SpoPlu   SpoRuf   SpoSch   SpoSim   SpoTor   SteDia   TacCor
## 1 3.483946 5.242597 4.825178 5.069198 5.544849 3.914820 41.60987 4.836622
## 2 3.227615 5.732869 4.247578 5.417510 5.466721 4.728831 41.75434 5.258086
## 3 3.497984 6.309418 2.992063 5.246964 5.149845 4.435311 42.64125 5.863616
## 4 3.149917 6.454133 2.951133 5.253867 5.456667 3.899500 44.07067 6.663067
## 5 3.446389 6.809431 3.171059 5.042793 5.492134 3.539630 45.07088 6.291722
## 6 3.981626 6.210236 3.725485 4.817861 5.919448 2.797162 46.70900 6.223065
##     TacCri   TacDel   TacLuc    TacPho   TacRuf   TacRuv   TacSur   TanArg
## 1 1.640398 8.762058 3.119652 0.7513315 8.583155 2.088650 1.789730 2.427432
## 2 1.815859 8.307054 3.777277 1.2612685 7.230387 1.740654 1.856085 2.181986
## 3 1.351065 7.036616 4.345204 1.2830036 4.991850 2.443249 1.607371 3.051023
## 4 1.128133 7.506867 3.852400 0.8894667 3.590933 3.159267 1.926733 2.819000
## 5 1.416974 7.649233 3.349419 0.5677444 2.922400 3.260626 1.982904 3.043611
## 6 1.414360 7.682982 3.427677 0.4146937 3.382674 3.840914 1.735825 3.130450
##     TanArt   TanCal   TanCay   TanChi   TanChr   TanCuc   TanCya   TanCyc
## 1 1.467633 3.939865 2.766577 1.546768 1.677360 5.353996 12.85256 24.99214
## 2 1.070534 4.167627 3.157926 2.042294 1.524777 4.903018 12.69680 25.51809
## 3 1.602883 4.549110 3.320213 2.138168 1.414596 5.572593 13.72460 26.39848
## 4 1.291167 4.865400 3.444867 2.011733 1.588752 6.040467 14.50273 26.02467
## 5 1.602079 5.802404 3.449726 2.003015 1.498061 5.336267 14.47962 26.77763
## 6 1.876052 6.399480 3.999852 2.592816 1.696240 6.013441 15.54088 28.33923
##     TanCyo   TanCyp     TanCyv   TanDes    TanDow   TanFas   TanFlo   TanFuc
## 1 4.685744 4.894128 1.53088288 3.585494  8.093189 7.288550 2.364691 4.202508
## 2 3.982295 5.122881 0.39092793 3.767771  7.414797 8.486483 2.275387 4.517622
## 3 4.159786 5.076450 0.01966066 3.400597  8.021703 9.367165 2.907186 4.599910
## 4 3.679571 5.279167 0.00000000 2.275361  7.958667 9.500778 3.105667 4.445889
## 5 3.799271 5.238824 0.00000000 2.683393  9.429824 8.347049 2.306123 5.335204
## 6 4.048019 5.988795 0.00000000 2.766825 11.963306 7.516703 2.475210 5.514745
##     TanGut   TanGyr   TanHei   TanIct   TanIno   TanJoh    TanLab   TanLar
## 1 4.829768 3.375108 4.754764 1.661829 4.501844 5.236825  8.742186 7.501444
## 2 4.057305 3.741441 4.464932 1.760623 4.582871 5.657435 10.544344 7.244984
## 3 3.251351 3.843132 3.552948 2.138634 4.196477 6.356735  9.801495 7.059277
## 4 3.480000 4.533933 3.714467 2.309267 4.342778 5.917067 10.277800 7.178000
## 5 3.814641 4.747278 1.882696 2.231774 4.564636 5.688391  9.468952 7.093171
## 6 3.383431 4.935061 1.599903 2.125413 4.895985 6.643610  9.604049 7.036554
##     TanLav   TanMex   TanMey   TanNic   TanNiv   TanPal   TanPar   TanPer
## 1 2.290441 2.827105 2.478378 21.84156 8.460319 11.34139 4.762243 5.691432
## 2 2.309212 3.017802 2.591824 22.41399 8.460668 11.56592 4.895393 6.298126
## 3 2.189450 2.984490 3.343730 23.14435 8.719510 12.82801 5.003978 6.256604
## 4 1.683000 3.178200 3.253833 23.55653 8.946267 13.70775 5.038467 6.772000
## 5 1.570694 3.372519 3.030491 24.11456 8.996733 13.53459 5.232363 6.668056
## 6 2.021892 3.450267 3.097725 25.70736 9.022232 14.43345 5.463470 5.572568
##      TanPre   TanPun   TanRue   TanRuf   TanRuu   TanSch    TanSel   TanVar
## 1  9.944220 6.008441 5.002063 3.424358 6.909399 1.863059 10.139715 3.795550
## 2  9.816647 6.111986 5.386405 2.702604 6.704295 1.906650  9.674646 3.480396
## 3 10.643818 6.772058 4.987523 2.554797 6.862164 1.649207  9.437048 3.359012
## 4 10.704467 7.040267 4.996000 2.483500 6.317833 1.751800 10.165000 3.968000
## 5 10.522515 7.286067 5.009926 2.402514 6.948602 1.717885 10.715284 4.519877
## 6 10.092454 7.242719 5.891712 2.773025 6.433117 1.826332 11.083459 5.116051
##     TanVas   TanVel   TanVir    TanVit   TanXac   TanXag   TerVir   ThlFul
## 1 2.090986 3.547241 5.851671 0.0000000 3.319324 5.198502 22.83156 6.217622
## 2 1.932000 3.829730 5.830550 0.0000000 3.283024 5.814634 23.48568 5.804049
## 3 1.741870 3.966005 5.949300 0.3065514 3.626273 6.880012 23.79442 6.021079
## 4 2.395533 4.065520 6.169000 0.3940000 3.822667 7.266778 25.07140 6.482333
## 5 3.007759 4.212550 6.293051 0.7635556 3.874549 7.289660 25.85069 6.294207
## 6 3.461418 4.419631 6.829633 0.8666216 4.072607 7.429174 26.23603 6.344178
##     ThlIno   ThlOrn   ThlPec   ThlRuf   ThlSor   ThrAbb   ThrBon   ThrCya
## 1 4.137212 4.018599 4.602153 4.298849 2.225132 16.54398 17.06939 23.94336
## 2 4.004766 4.506939 5.335471 4.226583 2.187126 16.95915 17.37653 23.97398
## 3 4.393178 5.093394 5.356084 4.089910 2.317160 17.08833 18.18349 24.66325
## 4 3.961167 5.019833 5.259889 3.864750 2.380533 17.96227 19.00967 26.03373
## 5 3.180060 4.355625 5.754611 3.825542 2.680167 18.30083 19.33483 26.94218
## 6 3.277047 4.366608 4.740420 3.410989 3.460623 18.89690 19.63418 27.52753
##     ThrCyp   ThrEpi    ThrGla   ThrOrn   ThrPal   ThrSay    TiaBic   TiaCan
## 1 11.77104 15.70287 11.164295 26.62071 11.45043 7.734886 0.8755788 5.133854
## 2 12.07851 15.92114 12.126169 26.56315 11.98017 7.490996 1.3769189 5.877162
## 3 11.92041 16.61791 11.888736 27.70014 12.96002 8.468750 2.1457500 6.069865
## 4 10.88247 17.18853 10.780833 29.30667 13.44800 9.263733 3.4594167 5.848733
## 5 11.06614 17.83998  9.608986 28.68717 13.98219 9.005744 3.1475093 5.876067
## 6 10.79618 17.93138  8.779264 28.67975 14.47977 9.377237 2.8914527 5.991153
##     TiaFul   TiaObs   TiaOli   TriMel   UroSto   VolJac   WetSte   XenCon
## 1 1.579715 4.562836 2.961937 4.663139 4.309291 4.797721 2.685617 6.996423
## 2 2.509440 4.996505 2.660923 4.703391 4.823276 4.793861 2.610338 7.475644
## 3 3.427998 6.191164 3.325404 5.140431 4.588535 5.222081 3.710387 6.674491
## 4 4.815200 5.576417 3.770933 5.100200 4.847111 5.361533 4.143833 6.251667
## 5 4.468637 6.331551 3.296244 4.975270 5.256043 5.136441 3.369602 5.716065
## 6 4.154216 6.832324 3.227438 4.348146 5.330523 4.315382 3.279829 5.451829
##     XenPar
## 1 11.79418
## 2 12.58418
## 3 12.37302
## 4 12.37817
## 5 13.14826
## 6 13.42591

Step 1 Extracting tanager tree

get the tanagers only

treepath <- "/Users/michaelalfaro/Dropbox/color_evolution_scratch/tanager_comp_analysis/allison_data/MCC_Tree_SpNames.nex"
tt <- read.nexus(treepath)

# Extract the tip labels (species names) from the tree
tan_species <- grep("^Tan", tt$tip.label, value = TRUE)

# Find the MRCA node of the species that start with "Tan"
mrca_node <- getMRCA(tt, tan_species)

# Extract the subtree from the MRCA node
tanagerTree <- extract.clade(tt, mrca_node)

# Plot the tanagerTree
p <- ggtree(tanagerTree) +
  geom_tiplab(size = 2.5, align = TRUE) +
  theme_tree2()

# Display the plot
print(p)


### Step 2 Create a Data Frame Matching the Tanager Tree


```r
# Extract the species names (tip labels) from the tanagerTree
tanager_species <- tanagerTree$tip.label

# Filter crown_data to include only species that are present in tanagerTree
tanData_filtered <- crown_data %>%
  select(wl, any_of(tanager_species))

# Display the filtered tanData to verify
head(tanData_filtered)
##    wl   TanCyo    TanLab   TanRue   TanChi   TanVel   TanCal   TanMex   TanIno
## 1 300 4.685744  8.742186 5.002063 1.546768 3.547241 3.939865 2.827105 4.501844
## 2 301 3.982295 10.544344 5.386405 2.042294 3.829730 4.167627 3.017802 4.582871
## 3 302 4.159786  9.801495 4.987523 2.138168 3.966005 4.549110 2.984490 4.196477
## 4 303 3.679571 10.277800 4.996000 2.011733 4.065520 4.865400 3.178200 4.342778
## 5 304 3.799271  9.468952 5.009926 2.003015 4.212550 5.802404 3.372519 4.564636
## 6 305 4.048019  9.604049 5.891712 2.592816 4.419631 6.399480 3.450267 4.895985
##      TanSel   TanFas   TanCyc     TanCyv   TanDes   TanChr   TanXac   TanSch
## 1 10.139715 7.288550 24.99214 1.53088288 3.585494 1.677360 3.319324 1.863059
## 2  9.674646 8.486483 25.51809 0.39092793 3.767771 1.524777 3.283024 1.906650
## 3  9.437048 9.367165 26.39848 0.01966066 3.400597 1.414596 3.626273 1.649207
## 4 10.165000 9.500778 26.02467 0.00000000 2.275361 1.588752 3.822667 1.751800
## 5 10.715284 8.347049 26.77763 0.00000000 2.683393 1.498061 3.874549 1.717885
## 6 11.083459 7.516703 28.33923 0.00000000 2.766825 1.696240 4.072607 1.826332
##     TanJoh   TanArt   TanIct   TanFlo   TanPar   TanLav   TanGyr   TanVas
## 1 5.236825 1.467633 1.661829 2.364691 4.762243 2.290441 3.375108 2.090986
## 2 5.657435 1.070534 1.760623 2.275387 4.895393 2.309212 3.741441 1.932000
## 3 6.356735 1.602883 2.138634 2.907186 5.003978 2.189450 3.843132 1.741870
## 4 5.917067 1.291167 2.309267 3.105667 5.038467 1.683000 4.533933 2.395533
## 5 5.688391 1.602079 2.231774 2.306123 5.232363 1.570694 4.747278 3.007759
## 6 6.643610 1.876052 2.125413 2.475210 5.463470 2.021892 4.935061 3.461418
##     TanNiv   TanFuc    TanDow   TanRuu   TanPun   TanXag   TanGut   TanVar
## 1 8.460319 4.202508  8.093189 6.909399 6.008441 5.198502 4.829768 3.795550
## 2 8.460668 4.517622  7.414797 6.704295 6.111986 5.814634 4.057305 3.480396
## 3 8.719510 4.599910  8.021703 6.862164 6.772058 6.880012 3.251351 3.359012
## 4 8.946267 4.445889  7.958667 6.317833 7.040267 7.266778 3.480000 3.968000
## 5 8.996733 5.335204  9.429824 6.948602 7.286067 7.289660 3.814641 4.519877
## 6 9.022232 5.514745 11.963306 6.433117 7.242719 7.429174 3.383431 5.116051
##     TanPal    TanVit   TanCuc   TanCay    TanPre   TanNic   TanCya   TanLar
## 1 11.34139 0.0000000 5.353996 2.766577  9.944220 21.84156 12.85256 7.501444
## 2 11.56592 0.0000000 4.903018 3.157926  9.816647 22.41399 12.69680 7.244984
## 3 12.82801 0.3065514 5.572593 3.320213 10.643818 23.14435 13.72460 7.059277
## 4 13.70775 0.3940000 6.040467 3.444867 10.704467 23.55653 14.50273 7.178000
## 5 13.53459 0.7635556 5.336267 3.449726 10.522515 24.11456 14.47962 7.093171
## 6 14.43345 0.8666216 6.013441 3.999852 10.092454 25.70736 15.54088 7.036554
##     TanCyp   TanVir   TanHei   TanArg   ThrCyp   ThrAbb   ThrPal   ThrOrn
## 1 4.894128 5.851671 4.754764 2.427432 11.77104 16.54398 11.45043 26.62071
## 2 5.122881 5.830550 4.464932 2.181986 12.07851 16.95915 11.98017 26.56315
## 3 5.076450 5.949300 3.552948 3.051023 11.92041 17.08833 12.96002 27.70014
## 4 5.279167 6.169000 3.714467 2.819000 10.88247 17.96227 13.44800 29.30667
## 5 5.238824 6.293051 1.882696 3.043611 11.06614 18.30083 13.98219 28.68717
## 6 5.988795 6.829633 1.599903 3.130450 10.79618 18.89690 14.47977 28.67975
##     ThrSay   ThrEpi   TanRuf
## 1 7.734886 15.70287 3.424358
## 2 7.490996 15.92114 2.702604
## 3 8.468750 16.61791 2.554797
## 4 9.263733 17.18853 2.483500
## 5 9.005744 17.83998 2.402514
## 6 9.377237 17.93138 2.773025

Step 3: Fitting Natural and Cubic Splines to the Tanager Data

### Step 3: Fitting Natural and Cubic Splines to the Tanager Data

# Function to fit a natural spline and a cubic spline, returning the models, coefficients, goodness-of-fit metrics, and RGB colors
fit_splines <- function(wl, reflectance) {
  # Fit natural cubic spline
  natural_spline_fit <- lm(reflectance ~ ns(wl, df = 5))
  
  # Fit cubic spline
  cubic_spline_fit <- lm(reflectance ~ poly(wl, 3, raw = TRUE))
  
  # Convert fitted splines to rspec and then to RGB hex codes
  natural_spline_rspec <- as.rspec(data.frame(wl = wl, reflectance = predict(natural_spline_fit)), lim = c(300, 700))
  natural_rgb <- spec2rgb(natural_spline_rspec)
  
  cubic_spline_rspec <- as.rspec(data.frame(wl = wl, reflectance = predict(cubic_spline_fit)), lim = c(300, 700))
  cubic_rgb <- spec2rgb(cubic_spline_rspec)
  
  # Extract coefficients
  natural_coef <- round(coef(natural_spline_fit), 2)
  cubic_coef <- round(coef(cubic_spline_fit), 2)
  
  # Calculate goodness-of-fit metrics
  natural_adj_r2 <- summary(natural_spline_fit)$adj.r.squared
  cubic_adj_r2 <- summary(cubic_spline_fit)$adj.r.squared
  
  natural_aic <- AIC(natural_spline_fit)
  cubic_aic <- AIC(cubic_spline_fit)
  
  return(list(
    natural_spline_fit = natural_spline_fit,
    cubic_spline_fit = cubic_spline_fit,
    natural_rgb = natural_rgb,
    cubic_rgb = cubic_rgb,
    natural_coef = natural_coef,
    cubic_coef = cubic_coef,
    natural_adj_r2 = natural_adj_r2,
    cubic_adj_r2 = cubic_adj_r2,
    natural_aic = natural_aic,
    cubic_aic = cubic_aic
  ))
}

# Apply the spline fitting function to each species
spline_results <- lapply(tanData_filtered[-1], function(reflectance) {
  fit_splines(tanData_filtered$wl, reflectance)
})
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 82 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## The spectral data contain 89 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 15 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 14 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 37 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 8 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 21 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 25 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 59 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 77 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 27 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 48 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 63 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.

## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## The spectral data contain 78 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 35 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 22 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 1 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 30 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 3 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 5 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1

Step 3b–making the species plots of natural and cubic spline

plot_splines_with_metrics <- function(species_name, wl, reflectance, spline_results, rgb_colors) {
  result <- spline_results[[species_name]]
  
  par(mfrow = c(1, 2))
  
  # Natural Spline
plot(wl, reflectance, type = "p", pch = 16, cex = 0.5,  # Adjusted size of points
     col = adjustcolor("black", alpha.f = 0.5),  # Adjusted transparency of points
     main = paste("Natural Spline -", species_name),
     xlab = "Wavelength (nm)", ylab = "Reflectance")
lines(wl, predict(result$natural_spline_fit), col = rgb_colors[[species_name]], lwd = 2)
mtext(paste("Adj R2:", round(result$natural_adj_r2, 3), "AIC:", round(result$natural_aic, 2)),
      side = 3, line = 0.5, cex = 0.8, col = "blue")

# Cubic Spline
plot(wl, reflectance, type = "p", pch = 16, cex = 0.5,  # Adjusted size of points
     col = adjustcolor("black", alpha.f = 0.5),  # Adjusted transparency of points
     main = paste("Cubic Spline -", species_name),
     xlab = "Wavelength (nm)", ylab = "Reflectance")
lines(wl, predict(result$cubic_spline_fit), col = rgb_colors[[species_name]], lwd = 2)
mtext(paste("Adj R2:", round(result$cubic_adj_r2, 3), "AIC:", round(result$cubic_aic, 2)),
      side = 3, line = 0.5, cex = 0.8, col = "blue")
}

for (species_name in names(spline_results)) {
  plot_splines_with_metrics(species_name, tanData_filtered$wl, tanData_filtered[[species_name]], spline_results, rgb_colors)
}

pdf("natural_vs_cubic_splines_with_metrics.pdf", width = 12, height = 8)
for (species_name in names(spline_results)) {
  plot_splines_with_metrics(species_name, tanData_filtered$wl, tanData_filtered[[species_name]], spline_results, rgb_colors)
}
dev.off()
## quartz_off_screen 
##                 2

Step 4: Plotting Spline Fits and Original Data in the Markdown Document

### Step 4: Calculate Average Splines and Standard Deviation

# Calculate the maximum reflectance across all splines
max_reflectance <- max(sapply(spline_results, function(result) {
  max(predict(result$natural_spline_fit), predict(result$cubic_spline_fit))
}))

# Initialize vectors to store average splines and standard deviations
natural_spline_avg <- numeric(length(tanData_filtered$wl))
cubic_spline_avg <- numeric(length(tanData_filtered$wl))
natural_spline_sd <- numeric(length(tanData_filtered$wl))
cubic_spline_sd <- numeric(length(tanData_filtered$wl))

# Calculate averages
for (result in spline_results) {
  natural_pred <- predict(result$natural_spline_fit)
  cubic_pred <- predict(result$cubic_spline_fit)
  
  natural_spline_avg <- natural_spline_avg + natural_pred
  cubic_spline_avg <- cubic_spline_avg + cubic_pred
}

# Convert sums to averages
natural_spline_avg <- natural_spline_avg / length(spline_results)
cubic_spline_avg <- cubic_spline_avg / length(spline_results)

# Calculate standard deviations
for (result in spline_results) {
  natural_pred <- predict(result$natural_spline_fit)
  cubic_pred <- predict(result$cubic_spline_fit)
  
  natural_spline_sd <- natural_spline_sd + (natural_pred - natural_spline_avg)^2
  cubic_spline_sd <- cubic_spline_sd + (cubic_pred - cubic_spline_avg)^2
}

natural_spline_sd <- sqrt(natural_spline_sd / (length(spline_results) - 1))
cubic_spline_sd <- sqrt(cubic_spline_sd / (length(spline_results) - 1))

### Correct for negative values in the average splines

# Convert the average natural and cubic splines to rspec objects
natural_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = natural_spline_avg), lim = c(300, 700))
## wavelengths found in column 1
cubic_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = cubic_spline_avg), lim = c(300, 700))
## wavelengths found in column 1
# Correct any negative values in the rspec objects
natural_spline_rspec <- procspec(natural_spline_rspec, fixneg = "zero")
## processing options applied:
## Negative value correction: converted negative values to zero
cubic_spline_rspec <- procspec(cubic_spline_rspec, fixneg = "zero")
## processing options applied:
## Negative value correction: converted negative values to zero
# Convert corrected rspec objects to RGB hex codes
natural_spline_rgb <- spec2rgb(natural_spline_rspec)
cubic_spline_rgb <- spec2rgb(cubic_spline_rspec)

### Step 4a: Plot Average Splines with Confidence Intervals

# Create a plot for the average natural spline with confidence intervals
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
     main = "Average Natural Spline with Variation Hotspots",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)), 
        c(natural_spline_avg + natural_spline_sd, rev(natural_spline_avg - natural_spline_sd)), 
        col = adjustcolor(natural_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)

# Create a plot for the average cubic spline with confidence intervals
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
     main = "Average Cubic Spline with Variation Hotspots",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)), 
        c(cubic_spline_avg + cubic_spline_sd, rev(cubic_spline_avg - cubic_spline_sd)), 
        col = adjustcolor(cubic_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)

# Save the hotspot plots to a PDF
pdf("spline_hotspots.pdf", width = 12, height = 8)
par(mfrow = c(1, 2))  # Plot both splines side by side

# Average Natural Spline with Hotspots in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
     main = "Average Natural Spline with Variation Hotspots",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)), 
        c(natural_spline_avg + natural_spline_sd, rev(natural_spline_avg - natural_spline_sd)), 
        col = adjustcolor(natural_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)

# Average Cubic Spline with Hotspots in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
     main = "Average Cubic Spline with Variation Hotspots",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)), 
        c(cubic_spline_avg + cubic_spline_sd, rev(cubic_spline_avg - cubic_spline_sd)), 
        col = adjustcolor(cubic_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)

dev.off()
## quartz_off_screen 
##                 2
### Step 4b: Plot Natural Spline Mean and Cubic Spline Mean with Species Overlaid

# Plot Natural Spline Mean with all species overlaid
par(mfrow = c(1, 2))

# Plot Natural Spline Mean with Species Overlaid
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
     main = "Natural Spline Mean with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
  species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$natural_spline_fit)), lim = c(300, 700))
  species_rspec <- procspec(species_rspec, fixneg = "zero")  # Correct any negative values
  species_rgb <- spec2rgb(species_rspec)
  lines(tanData_filtered$wl, predict(spline_results[[i]]$natural_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 82 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 37 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 25 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 27 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 63 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 30 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)

# Plot Cubic Spline Mean with Species Overlaid
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
     main = "Cubic Spline Mean with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
  species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$cubic_spline_fit)), lim = c(300, 700))
  species_rspec <- procspec(species_rspec, fixneg = "zero")  # Correct any negative values
  species_rgb <- spec2rgb(species_rspec)
  lines(tanData_filtered$wl, predict(spline_results[[i]]$cubic_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 89 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 15 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 14 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 8 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 21 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 59 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 77 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 48 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 78 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 35 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 22 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 1 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 3 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 5 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)

# Save the overlay plots to a PDF
pdf("spline_mean_with_species_overlay.pdf", width = 12, height = 8)
par(mfrow = c(1, 2))

# Plot Natural Spline Mean with Species Overlaid in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
     main = "Natural Spline Mean with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
  species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$natural_spline_fit)), lim = c(300, 700))
  species_rspec <- procspec(species_rspec, fixneg = "zero")  # Correct any negative values
  species_rgb <- spec2rgb(species_rspec)
  lines(tanData_filtered$wl, predict(spline_results[[i]]$natural_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 82 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 37 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 25 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 27 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 63 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 30 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)

# Plot Cubic Spline Mean with Species Overlaid in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
     main = "Cubic Spline Mean with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
  species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$cubic_spline_fit)), lim = c(300, 700))
  species_rspec <- procspec(species_rspec, fixneg = "zero")  # Correct any negative values
  species_rgb <- spec2rgb(species_rspec)
  lines(tanData_filtered$wl, predict(spline_results[[i]]$cubic_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 89 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 15 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 14 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 8 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 21 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 59 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 77 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 48 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 78 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 35 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 22 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 1 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 3 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## The spectral data contain 5 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)

dev.off()
## quartz_off_screen 
##                 2

Step 4b Make a pair of plots of the natural and cubic spline where the spline is colored as a spec that is converted to its hsv value (make sure no negative values)

### Step 4b: Plot Natural Spline with Species Overlaid and Print to PDF

# Convert the natural spline to an rspec object and then to an RGB hex code
natural_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = natural_spline_avg), lim = c(300, 700))
## wavelengths found in column 1
natural_spline_rgb <- spec2rgb(natural_spline_rspec)

# Convert the cubic spline to an rspec object and then to an RGB hex code
cubic_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = cubic_spline_avg), lim = c(300, 700))
## wavelengths found in column 1
cubic_spline_rgb <- spec2rgb(cubic_spline_rspec)

# Plot the natural and cubic splines side by side in R Markdown
par(mfrow = c(1, 2))  # Arrange plots side by side

# Plot Natural Spline with Species Overlaid
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,  # Increased thickness
     main = "Natural Spline with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))

# Loop through each species to plot their splines with slightly more visible transparency
for (species in names(spline_results)) {
  species_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$natural_spline_fit)), lim = c(300, 700))
  species_rgb <- spec2rgb(species_spline_rspec)
  
  # Adjust transparency to make species lines slightly more visible
  species_rgb_alpha <- adjustcolor(species_rgb, alpha.f = 0.4)  # Adjusted alpha for more visibility
  
  lines(tanData_filtered$wl, predict(spline_results[[species]]$natural_spline_fit), col = species_rgb_alpha)
}
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 82 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 37 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 25 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 27 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 63 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.

## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 30 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
# Re-plot the natural spline on top to ensure it stands out
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)  # Increased thickness

# Plot Cubic Spline with Species Overlaid
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,  # Increased thickness
     main = "Cubic Spline with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))

# Loop through each species to plot their splines with slightly more visible transparency
for (species in names(spline_results)) {
  species_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$cubic_spline_fit)), lim = c(300, 700))
  species_rgb <- spec2rgb(species_spline_rspec)
  
  # Adjust transparency to make species lines slightly more visible
  species_rgb_alpha <- adjustcolor(species_rgb, alpha.f = 0.4)  # Adjusted alpha for more visibility
  
  lines(tanData_filtered$wl, predict(spline_results[[species]]$cubic_spline_fit), col = species_rgb_alpha)
}
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 89 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 15 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 14 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 8 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 21 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 59 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 77 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 48 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 78 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 35 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 22 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 1 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 3 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 5 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
# Re-plot the cubic spline on top to ensure it stands out
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)  # Increased thickness

# Save the plots to a PDF as well
pdf("natural_and_cubic_splines_with_species_overlay.pdf", width = 12, height = 8)
par(mfrow = c(1, 2))  # Arrange plots side by side in the PDF

# Plot Natural Spline with Species Overlaid in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
     main = "Natural Spline with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (species in names(spline_results)) {
  species_rgb_alpha <- adjustcolor(spec2rgb(as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$natural_spline_fit)), lim = c(300, 700))), alpha.f = 0.4)
  lines(tanData_filtered$wl, predict(spline_results[[species]]$natural_spline_fit), col = species_rgb_alpha)
}
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 82 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.

## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 37 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 25 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 27 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 63 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.

## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 30 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)

# Plot Cubic Spline with Species Overlaid in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
     main = "Cubic Spline with Species Overlaid",
     xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (species in names(spline_results)) {
  species_rgb_alpha <- adjustcolor(spec2rgb(as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$cubic_spline_fit)), lim = c(300, 700))), alpha.f = 0.4)
  lines(tanData_filtered$wl, predict(spline_results[[species]]$cubic_spline_fit), col = species_rgb_alpha)
}
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 89 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 15 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 14 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 8 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 21 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 59 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 77 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 48 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 78 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 35 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 22 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## The spectral data contain 1 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 3 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## The spectral data contain 5 negative value(s),
##  which may produce unexpected results if used in models.
## Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)

dev.off()
## quartz_off_screen 
##                 2

Step 5 Print a Summary Table of Spline Fits

# Open the PDF device
pdf("spline_fits_with_rgb_and_metrics.pdf", width = 12, height = 8)

# Set up the plotting area to show two plots per page
par(mfrow = c(1, 2))

# Loop over each species and plot the original data and spline fits side by side
for (species in names(spline_results)) {
  result <- spline_results[[species]]
  
  # Get the correct RGB color
  color <- rgb_colors[grep(species, names(rgb_colors))]
  
  # Plot original data and natural spline fit
  plot(tanData_filtered$wl, tanData_filtered[[species]], main = paste("Natural Spline Fit for", species),
       xlab = "Wavelength (nm)", ylab = "Reflectance", type = "o", col = color)
  lines(tanData_filtered$wl, predict(result$natural_spline_fit), col = "red", lwd = 2)
  mtext(paste("Adj R2:", round(result$natural_adj_r2, 3), "AIC:", round(result$natural_aic, 2)),
        side = 3, line = 0.5, cex = 0.7, col = "blue")
  
  # Plot original data and cubic spline fit
  plot(tanData_filtered$wl, tanData_filtered[[species]], main = paste("Cubic Spline Fit for", species),
       xlab = "Wavelength (nm)", ylab = "Reflectance", type = "o", col = color)
  lines(tanData_filtered$wl, predict(result$cubic_spline_fit), col = "green", lwd = 2)
  mtext(paste("Adj R2:", round(result$cubic_adj_r2, 3), "AIC:", round(result$cubic_aic, 2)),
        side = 3, line = 0.5, cex = 0.7, col = "blue")
}

# Close the PDF device
dev.off()
## quartz_off_screen 
##                 2

Step 6 Print a Summary Table of Spline Fits

# Create the spline summary table with coefficients and goodness-of-fit metrics
spline_table <- data.frame(
  Species = names(spline_results),
  Natural_Spline_Coefficients = sapply(spline_results, function(x) paste(x$natural_coef, collapse = ", ")),
  Cubic_Spline_Coefficients = sapply(spline_results, function(x) paste(x$cubic_coef, collapse = ", ")),
  Natural_Adj_R2 = sapply(spline_results, function(x) round(x$natural_adj_r2, 3)),
  Cubic_Adj_R2 = sapply(spline_results, function(x) round(x$cubic_adj_r2, 3)),
  Natural_AIC = sapply(spline_results, function(x) round(x$natural_aic, 2)),
  Cubic_AIC = sapply(spline_results, function(x) round(x$cubic_aic, 2))
)

# Remove row names
rownames(spline_table) <- NULL

# Display the table in the document without row names
kable(spline_table, caption = "Summary of Spline Fits for Each Species (Including Goodness-of-Fit)") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Summary of Spline Fits for Each Species (Including Goodness-of-Fit)
Species Natural_Spline_Coefficients Cubic_Spline_Coefficients Natural_Adj_R2 Cubic_Adj_R2 Natural_AIC Cubic_AIC
TanCyo 5.49, 1.85, 6.47, 0.09, 2.03, 0.44 9, -0.07, 0, 0 0.920 0.717 434.02 940.85
TanLab 8.82, -2.1, 8.11, 1.57, -3.02, 5.53 124.37, -0.81, 0, 0 0.971 0.769 621.19 1447.47
TanRue 3.55, 0.23, 29.51, 5.17, 10.86, 2.65 177.07, -1.28, 0, 0 0.970 0.770 1229.01 2045.23
TanChi -4.73, -26.22, 79.79, 28.04, 39.68, -0.95 902.87, -6.16, 0.01, 0 0.867 0.659 2779.55 3156.06
TanVel 4.67, 3.79, -0.58, 0.53, -2.32, 2.04 -25.97, 0.17, 0, 0 0.854 0.317 419.45 1034.89
TanCal 5.14, 20.49, 66.62, 86.72, 56.18, 37.57 812.16, -5.72, 0.01, 0 0.995 0.984 1739.27 2159.82
TanMex 2.83, 1.49, 1.72, 1.33, 3.24, 1.8 -13.24, 0.1, 0, 0 0.851 0.848 -34.11 -29.58
TanIno 5.31, 4.18, 4.8, 5.56, 8.88, 6.11 -17.72, 0.13, 0, 0 0.983 0.983 82.00 89.03
TanSel 2.54, -4.87, 117.83, 13.75, 31.09, -7.3 783.68, -5.71, 0.01, 0 0.925 0.715 2730.48 3265.24
TanFas 1.31, 2.88, 112.54, 10.91, 28.07, 0.89 634.65, -4.76, 0.01, 0 0.956 0.743 2471.70 3176.45
TanCyc 37.56, 12.07, -18.39, -27.75, -32.97, -24.42 -323.85, 2.31, 0, 0 0.969 0.911 1788.56 2206.97
TanCyv 0.77, -4.97, 21.78, 29.97, 37.82, 27.44 274.9, -1.78, 0, 0 0.985 0.943 1483.02 2024.56
TanDes -1.07, -19.37, 66.97, 2.65, 28.12, -0.49 509.47, -3.53, 0.01, 0 0.902 0.550 2331.92 2942.65
TanChr 1.29, 3.58, 3.84, 8.51, 6.95, 3.27 34.67, -0.25, 0, 0 0.978 0.931 260.33 713.45
TanXac 3.57, -0.7, 2.63, 59.27, 56.46, 50.7 301.81, -1.81, 0, 0 0.995 0.939 1482.10 2490.76
TanSch -0.57, -9.88, 39.78, 44.33, 57.75, 42.84 471.6, -3.1, 0.01, 0 0.980 0.935 1968.85 2440.76
TanJoh 1.93, -14.12, 49.87, 12.8, 25.72, 3.14 474.41, -3.23, 0.01, 0 0.920 0.678 2112.46 2668.77
TanArt 0.23, -3.98, 13.2, 36.54, 41.89, 35.77 250.37, -1.57, 0, 0 0.993 0.960 1345.99 2017.60
TanIct -0.48, -9.35, 34, 33.51, 43.72, 27.5 410.4, -2.71, 0.01, 0 0.971 0.904 1889.41 2366.07
TanFlo -1.53, -24.02, 70.46, 59.26, 72.43, 41.1 927.96, -6.16, 0.01, 0 0.960 0.871 2512.40 2981.26
TanPar 3.34, -7.99, 33.57, 53.71, 53.47, 42.64 559.29, -3.66, 0.01, 0 0.993 0.959 1595.16 2304.96
TanLav 2.05, 0.16, -0.32, 12.54, 24.61, 32.5 -58.09, 0.49, 0, 0 0.999 0.996 50.66 697.53
TanGyr 4.04, 0.37, -0.54, 10.6, 19.7, 25.76 -40.58, 0.37, 0, 0 0.997 0.992 416.94 824.40
TanVas 4.47, 27.81, 20.65, 17.18, 10.92, 22.3 -98.7, 0.4, 0, 0 0.979 0.817 1375.32 2239.01
TanNiv 6.08, 9.48, 108.33, 41.1, 29.4, 15.17 875.62, -6.37, 0.01, 0 0.989 0.883 2025.91 2956.90
TanFuc 4.26, 3.87, 44.72, 14.14, 15.7, 12.45 297.85, -2.18, 0.01, 0 0.995 0.861 939.51 2282.45
TanDow 9.29, -4.15, 5.65, 8.11, 4.92, 14.15 149.15, -0.92, 0, 0 0.987 0.972 731.23 1026.97
TanRuu 5.09, 1.49, 3.27, 5.51, 6.82, 8.33 17.61, -0.09, 0, 0 0.948 0.947 739.21 743.41
TanPun 3.82, 7.76, 35.69, 28.25, 22.32, 5.18 321.87, -2.3, 0.01, 0 0.955 0.919 1806.05 2040.05
TanXag 5.57, -29.14, 57.48, 4.05, 24.95, -8.26 617.88, -4.1, 0.01, 0 0.941 0.518 2071.70 2908.52
TanGut 2.26, -8.05, 35.49, 12.61, 23.07, 3.8 319.67, -2.17, 0, 0 0.956 0.737 1607.22 2318.38
TanVar 1.57, -13.14, 39.81, 4.7, 21.3, -2.47 331.05, -2.24, 0, 0 0.940 0.561 1717.09 2512.59
TanPal 15.54, 3.54, 5.36, 0, 14.83, -2.49 -102.78, 0.78, 0, 0 0.863 0.605 891.24 1315.52
TanVit 2.25, 2.13, 5.97, 38.59, 39.52, 38.61 168.59, -1.03, 0, 0 0.999 0.969 661.58 1912.92
TanCuc 5.49, 0.1, 0.77, 2.67, 5.01, 6.22 -0.38, 0.05, 0, 0 0.990 0.990 -245.67 -237.98
TanCay 3.33, 0.12, 8.37, 27.98, 28.83, 29.58 165.23, -1.03, 0, 0 0.999 0.984 107.89 1436.42
TanPre 9.87, -6.53, 1.96, 43.72, 25.99, 27.87 430.45, -2.69, 0.01, 0 0.995 0.938 1169.47 2211.25
TanNic 26.26, 42.67, 23.41, 0.32, 46.88, -15.12 -680.49, 4.43, -0.01, 0 0.994 0.989 1051.49 1327.63
TanCya 10.54, 0.71, 81.85, -10.46, 53.36, -36.27 47.12, -0.51, 0, 0 0.977 0.661 1882.47 2963.59
TanLar 9.69, 14.72, 14.22, 74.71, 67.17, 80.11 288.47, -1.82, 0, 0 0.989 0.975 2016.22 2358.96
TanCyp 5, 0.16, 0.48, 1.2, 1.1, 2.38 10.89, -0.04, 0, 0 0.904 0.893 -49.47 -6.38
TanVir 6.5, -0.97, -0.72, 0.13, -1.12, 1.66 22.06, -0.1, 0, 0 0.868 0.824 -105.50 7.12
TanHei 3.85, 0.96, 2.17, 1.88, 3.42, 2.48 0.76, 0.01, 0, 0 0.767 0.762 500.96 507.78
TanArg 2.59, 0.24, 0.56, 0.82, 0.87, 0.74 6.28, -0.03, 0, 0 0.616 0.616 -28.33 -30.27
ThrCyp 10.82, -3.61, 17.18, -0.92, 9.05, -2.02 92.98, -0.58, 0, 0 0.955 0.485 812.06 1786.58
ThrAbb 11.25, 7.92, 4.52, -11.19, 46.2, -25.53 -503.23, 3.51, -0.01, 0 0.980 0.780 1535.84 2493.80
ThrPal 12.04, -16.38, 18.35, -8.19, 43.63, -15.23 -69.59, 0.72, 0, 0 0.950 0.238 1581.24 2670.14
ThrOrn 26.98, 1.08, -10.46, -15.61, 4.52, -21.38 -258.72, 1.97, 0, 0 0.982 0.930 1290.36 1821.61
ThrSay 7.51, 5.2, 18.59, 11.41, 24.03, 1.39 28.85, -0.2, 0, 0 0.976 0.783 732.42 1614.16
ThrEpi 17.47, 7.34, 36.49, -1.53, 31.31, -2.88 -83.1, 0.53, 0, 0 0.969 0.664 1193.19 2149.73
TanRuf 2.9, 2.46, 8.23, 15.08, 17.93, 19.02 63.51, -0.41, 0, 0 0.999 0.997 -100.03 243.82
# Save the table as a CSV
write.csv(spline_table, "spline_summary_table.csv", row.names = FALSE)

# Fixing the PDF export by increasing the margins and adjusting the layout
pdf("spline_summary_table.pdf", width = 11, height = 8.5)
gridExtra::grid.table(spline_table, rows = NULL)
dev.off()
## quartz_off_screen 
##                 2

Step 7 Save data objects for beast run

# Step 7: Save the Tanager Tree and Data as an R Data Object

# Save the tree, filtered spectral data, and spline results into an RData file
save(tanagerTree, tanData_filtered, spline_results, file = "tanager_analysis_results.RData")

# Inform the user that the data has been saved
cat("The Tanager tree, data, and spline results have been saved to 'tanager_analysis_results.RData'.\n")
## The Tanager tree, data, and spline results have been saved to 'tanager_analysis_results.RData'.

Step 8 Getting things ready for BEAST

# Export spline coefficients as a CSV file for BEAST
spline_traits <- data.frame(Species = names(spline_results))

# Add the coefficients as traits
for (i in 1:length(spline_results)) {
  spline_traits[i, "Natural_Spline_1"] = spline_results[[i]]$natural_coef[2]
  spline_traits[i, "Natural_Spline_2"] = spline_results[[i]]$natural_coef[3]
  spline_traits[i, "Natural_Spline_3"] = spline_results[[i]]$natural_coef[4]
  spline_traits[i, "Natural_Spline_4"] = spline_results[[i]]$natural_coef[5]
  spline_traits[i, "Cubic_Spline_1"] = spline_results[[i]]$cubic_coef[2]
  spline_traits[i, "Cubic_Spline_2"] = spline_results[[i]]$cubic_coef[3]
  spline_traits[i, "Cubic_Spline_3"] = spline_results[[i]]$cubic_coef[4]
}

# Save the spline traits to a CSV file
write.csv(spline_traits, "spline_coefficients_for_BEAST.csv", row.names = FALSE)

Step 9 BEAST xml

not sure yet how to do this but here is the general appropach….

Load Tree and Data: Load your phylogenetic tree (e.g., tanagerTree) and the CSV file containing the traits. Define the Traits: In the XML file, define each spline coefficient as a separate trait. If you want to estimate a single evolutionary rate for all coefficients, you can define them as independent traits. Model Specification: For each trait, specify a Brownian motion model (bm). If you want to treat the coefficients independently but estimate the same evolutionary rate for all, use a shared rate parameter across the traits. Rate Estimation: Define a hyperparameter for the rate of evolution and link it to all the traits. Alternatively, run separate BEAST analyses for each trait if you want to estimate the rates independently.

<trait id="natural_spline_1" spec="RealParameter" dimension="N">
    <statefile name="state" startstate="true"/>
</trait>

<distribution id="prior" spec="util.CompoundDistribution">
    <distribution spec="beast.math.distributions.MRCAPrior" id="TMRCA">
        <tree topologies="$tanagerTree" spec="TreeLikelihood">
            <distribution spec="beast.evolution.likelihood.MG94Likelihood" id="treeModel">
                <data id="natural_spline_1.data" spec="beast.evolution.alignment.Alignment" dataType="continuous">
                    <sequence taxon="species1" totalcount="4" value="coeff1, coeff2, coeff3, coeff4"/>
                    <sequence taxon="species2" totalcount="4" value="coeff1, coeff2, coeff3, coeff4"/>
                </data>
            </distribution>
        </tree>
    </distribution>
</distribution>

<run id="mcmc" spec="MCMC">
    <state>
        <parameter id="rate" name="stateNode" dimension="1">0.1</parameter>
    </state>

    <posterior>
        <distribution id="posterior" spec="util.CompoundDistribution">
            <prior id="prior" spec="beast.math.distributions.MRCAPrior" tree="@treeModel">
                <normalPrior mean="0" stdev="1" name="distr"/>
            </prior>
            <likelihood id="likelihood" spec="util.CompoundDistribution">
                <distribution id="treeLikelihood" spec="TreeLikelihood">
                    <data id="alignment" spec="Alignment" dataType="continuous"/>
                    <siteModel id="SiteModel" spec="SiteModel" gammaCategoryCount="4">
                        <substitutionModel id="F81" spec="F81"/>
                    </siteModel>
                </distribution>
            </likelihood>
        </distribution>
    </posterior>
</run>

You will need to replicate this setup for each coefficient, ensuring that they either share or have independent rate parameters depending on your choice.

  1. Run BEAST

Load the XML into BEAST: Once the XML is set up, load it into BEAST. Run the Analysis: Start the analysis to estimate the evolutionary rates. 4. Post-Analysis

After the BEAST run completes, you can analyze the results using Tracer to check the estimated rates and other parameters.